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Abstract 

We present numerical enhancements of a multiscale domain decomposition strategy based on 
a LaTIn solver and dedicated to the computation of the debounding in laminated composites. We 
^ show that the classical scale separation is irrelevant in the process zones, which results in a drop 

^ ^ in the convergence rate of the strategy. We show that performing nonlinear subresolutions in the 

vicinity of the front of the crack at each prediction stage of the iterative solver permits to restore 
the effectiveness of the method. 
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^ Introduction 

^ The reliable simulation of delamination in laminated composites requires fine models designed at the 

micro scale [T4| (where fibers can be distinguished) or at most at the meso scale [19l [3] (where plies 
^ can be distinguished), because at these levels physics can be correctly taken into account. Thus, 

I— I even the simulation of small laminated components implies to use huge discrete models which are out 

^ of the reach of non-optimized computational techniques. In this paper we mainly focus on two key 

^ ingredients to set up efficient strategies: the handling of nonlinearity and the use of nested scales of 

^\ calculation to quickly distribute the computational effort on the whole structure — or more precisely on 

On all substructures since using domain decomposition is almost mandatory to achieve high-performance 

C parallelism. 

Most classical strategies provide separated answers to these two problems: nonlinearity is handled 
On with Newton-Raphson algorithm (with, if required, arc-length control) [17| or with asymptotic numeri- 

cal method [24] ; multiscaling is realized on the linearized systems with domain decomposition methods 
^ ] like FETI [5], BDD [20], Schwarz [HI [El [2]. Unfortunately, it is well known that the convergence of 

. . such approaches can be seriously impaired in the case of strong localized non-linearities [1] : not only 

^ the nonlinear process may require lots of load increments to converge but linear systems may be poorly 

conditioned and difficult to solve. 
^ Recent studies have tried to take advantage of the domain decomposition to contain the difficulties 

^ associated to the nonlinearity within subdomains: they are based a process called nonlinear relocaliza- 

tion [5j which consists in solving nonlinear problems independently inside subdomains (with well-chosen 
interface conditions). In ^22J this procedure was interpreted as conducting a Newton-Raphson on a 
nonlinear condensed problem, which enabled to classify variants and to propose new algorithms. Stud- 
ies in [7] somehow consist in applying the relocalization philosophy under specific software constraints 
(use of a commercial "closed" finite element software). 

Our studies are based on the LaTIn method which, from the very beginning [llj , has been designed 
as a nonlinear solver where nonlinearity was dealt with at the smallest possible scale (typically inde- 
pendently on each Gauss point for local material constitutive laws). The method was then adapted to 
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substructuring by the introduction of unknown kinematic (displacement) and static (traction) interface 
fields which were linked by a constitutive equation (perfect joint, elastic joint, contact, friction). The 
Schwarz-type resolution algorithm was highly parallel, main operations were: resolution of nonlinear 
problem at the Gauss-point-scale, resolution of sparse linear systems at the subdomain-scale, exchange 
of interface vectors between subdomains [12,. Yet that local handling of information (subdomains 
communicating only with their neighbors) lead to a non-scalable algorithm: for a given problem, the 
convergence rate decreased when the number of subdomains increased. Scalability was achieved using 
a multiscale extension which consisted in insuring partial continuity and equilibrium conditions be- 
tween subdomains which resulted in the resolution of a small global (defined on the whole structure) 
linear problem 113j. The weak continuity conditions (called "macro" conditions) were inspired from 
Saint- Venant principle and homogenization techniques 23 , 6 : the long-range influence of the physical 
phenomena was thus transmitted to the whole structure, so that subdomains not only got information 
from their neighbors but also from distant substructures. The number of iterations to converge thus 
became independent on the number of subdomains. Since then the method has been validated on 
various nonlinear problems [211 US] • 

Because of the presence of both kinematic and static fields on the interfaces, introducing cohesive 
behaviors [Tj on the interfaces is very easy in the LaTIn method 10 . Unfortunately for such interface 
behaviors, insuring the scalability of the method is not as straightforward: the stress singularity 
associated to the crack is not well captured within classical macro quantities which results in the 
long range effect not being transmitted and the convergence being seriously impaired (see also |9j for 
a somehow similar study). One interpretation of this phenomenon is that the very local treatment 
of nonlinearity prevent us from filtering the long range information to be transmitted globally. In 
this paper, we show that an intermediate treatment of the nonlinearity enables to detect the relevant 
information to spread along the structure. This technique is thus connected to the previously mentioned 
relocalization techniques, though here the treatment of nonlinearity is scaled-up instead of being scaled- 
down which seems more robust with respect to the risk of non-physical local instabilities, moreover 
very pertinent boundary conditions are easily imposed on the boundary of the relocalization domain. 

The rest of the paper is organized as follow: in the first section, we present the classical LaTIn 
method; in the second section, we show that damaging cohesive behavior leads to the loss of the 
scalability of the method and that strategies based on the a priori enrichment of the macro information 
to transmit are doomed to inefficiency; the third section presents the relocalization technique; the 
fourth section opens the discussion on the method. 

1 Reference debounding problem and resolution strategy 
1.1 Substructuring of the debounding problem 

The laminated structure E occupying the domain fl is made out of adjacent plies, separated by cohesive 
interfaces. An external traction field (respectively displacement field U_^) is prescribed on Part dilf 
(respectively dflu) of the boundary dfl. The volume force is denoted The simulation is performed 
under the assumption of small perturbations and the evolution over time is supposed to be quasi-static 
and isothermal; classical incremental scheme is used. 

The structure is decomposed into substructures and interfaces as represented Fig. [l] Each of these 
mechanical entities has its own kinematic and static unknown fields as well as its own constitutive 
law. The substructuring is driven by the will to match domain decomposition interfaces with material 
cohesive interfaces, so that each substructure E belongs to a unique ply and has a constant linear 
constitutive law. Let be the Cauchy stress tensor in E, the displacement field and €{ue) the 
symmetric part of the displacement gradient. A substructure E defined in Domain £7^; is connected 
to an adjacent substructure E' through an interface Tee' — dflE H dil,E' (Fig-[2|. The surface entity 
Tee' applies force distributions Fe, Fe' 'well as displacement distributions W_e, W-E' to E and E' 
respectively. Let us denote Te ~ Ub'se -'^■EB' ■ On a substructure E such that O^Ie n d^l ^ 0, the 
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Perfect interfaces 
Figure 1: Substructuring of the composite structure 

boundary condition (C/^,_F^) is applied through a boundary interface Te^- 




Figure 2: Mixed description of the unknown fields 

At each step of the incremental time resolution scheme, the substructured quasi-static problem 
consists in finding s = {sE)EeE, where se — iW-E^E-E)^ which is a solution to the following equations: 



• Kinematic admissibility of Substructure E: 

at each point of F^ , Ue — We 

• Static admissibility of Substructure E: 

y{uE*,WE*) EUeXWe / MfiVn. = ^E*, 



Tr 



'-E 



f l^.UE*dn+ f FeWe^cIT 



• Linear orthotropic constitutive law of Substructure E: 

at each point of 51^;, = Ke(M^) 



(1) 



(2) 



(3) 
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• Behavior of each interface Tee'- 



at each point oiTEE', Hee' {Ke , , , Zb' ) = (4) 
• Behavior of the interfaces at the boundary dil.f n F^: 

at each point of Te^, TlEai^EiPE) = /r) 
(We = Ud on d^lu and Fe = Z^ '^^ 9^f) 

The formal relation Ree' = is now made explicit in two representative cases: 

Perfect interface: < 77^' ^ 4"^' (6) 
i BLe ~ B^E' =0 

f Fe + Fe,=0 

Cohesive interface: < ^ ^, fnjin ^ A mn (7) 

The last equation is the nonlinear constitutive law of the cohesive interfaces. The progressive softening 
of these interfaces (from healthy elastic to zero-stiffness in traction and shear) with respect to the 
history of the interface variables is described using continuum damage mechanics. Further details on 
the model used and on identification issues can be found in 

1.2 Two-scale iterative resolution of the substructured problem 

1.2.1 Introduction of the macroscopic scale 

The substructured problem defined in previous section shall eventually be solved by an iterative LaTIn 
algorithm, which will be detailed in the next subsection. In order the strategy to be scalable (that is, 
for a given problem, to converge independently on the substructuring), a global coarse grid problem 
must be solved at each iteration of the solver. This coarse problem is associated to the equilibrium 
and continuity of so-called "macroscopic" force and displacement fields of the interfaces. 

On each interface Tee' such that {E,E') S E^, the interface fields are split into a macro part 
and a micro part the former belonging to a small-dimension subspace (9 macro degrees of freedom 
per plane interface in 3D). 

Fe = F^' + F!^ We = W^ + W^ (8) 

The macro and micro data are uncoupled with respect to the interface virtual work: 

y{FE,WE) ^TexWe, 

Fe.WecIT^ f F%^.W^' dT+ f F'J^.W'J^ dV 

Macro spaces are determined by the choice of their basis. Numerical tests have shown that using 
a linear macro basis provides, in general, good scalability properties to the method [[ISj]. Indeed, 
the corresponding macro space includes the part of the interface fields with the highest wavelength. 
Consequently, according to the Saint- Venant principle, the micro complement found iteratively through 
the resolution of local problems only has a local influence. 

1.2.2 The iterative algorithm 

The iterative LaTIn algorithm [12;], designed to solve nonlinear problems, is here applied to the 
resolution of the substructured debounding problem, the nonlinearities being lumped in the (cohesive) 
interfaces. 
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The equations of the problem can be split into the set of linear equations in substructure and inter- 
face macroscopic variables (static and kinematic admissibility of the substructures, linear constitutive 
law of the substructures, linear equilibrium of the macro interface forces) and the set of local equations 
in interface variables (behavior of the interfaces). The solutions s {se)e<£'e = {]^e, £^e)e£E to the 
first set of equations belong to Space Ad, while the solutions s = i'sE)EeE — (W^^ , f ^ ) kpf, to the 
second set of equations belong to T. Hence, the converged solution Sref is such that Sref & A.df]T. 




Figure 3: Illustration of the LaTIn iterative algorithm 

The resolution scheme consists in searching for the solution Sref alternatively in these two spaces 
along search directions E+ and (see Fig. wl): 



• Find Sfj+i G T such that ySn+^ ^ ■^nj G E+ (local stage) 

• Find Sn+i G Ad such that (^s„+i — 1 ^ G (linear stage) 
In the following, the subscript n will be dropped. 

Local stage Independent local problems are solved at each point of the interfaces {TEE')\(E,E')eE^ 
(and {TEa)EeE for the interfaces belonging to 9r2„ U d^lf): 

nEE'{WE,WE,,FE,FE,)=0 

Find {Fe, We,Fe,, We,) such that: { [f^^Fj^)- k+iW^ - VKb) = (10) 

{lE'-FE,)-k+(WE,-WE,)^0 



The last two equations are the search direction equations E+. For a cohesive interface. Problem (10) 
is nonlinear, and solved by a Newton-Raphson scheme. 

Linear stage The linear stage consists in solving linear systems in substructure variables under the 
constraint of macroscopic equilibrium of the interface forces: 



V Tee' I {E, E') G E^, F*^ + F% = 

V r^,, / n a% ^ 0, z*^-zf = o 



(11) 



The macroscopic condition is not compatible with the monoscale search direction E~ coupling the 
interface displacement and forces fields at the linear stage. Hence the search direction is weakened and 
verified at best under the macroscopic constraint [|15)]. Technically this is realized using a Lagrangian 
whose stationarity leads to a modified local search direction: 

\/We* &We, I (FE-pE + k' {WE-^E)-k'W.^^) ■We'' dV = Q (12) 
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where the Lagrange multiplier is a macroscopic unknown of Interface T ee' 



The expression of the problem to solve on each substructure E is obtained by substituting ( 12 ) in 



"iiiiE^WE") e UexWe, 



l^.UE*dn+ / {F 



k- 

krw 



We -We* dT 



(13) 



-A-r 



k-W ).We* dV 



M 



The condensation of this equation on the macro degrees of freedom leads to a relation coupling 



and W_E which can be introduced in the macro equilibrium equation (11). Eventually, one gets a 
small linear system defined on macro degrees of freedom. All subdomains contribute to that "global" 
system through explicitly computed homogenized (condensed) flexibilities L^^. 



VW e , ^ / Lf .li^ dV 

E •^'"e 

= E / ^d-w dr-J2 Fe-W dv 

c; JTEndnf £, JTe 



(14) 



The right-hand side of Equation ( 14 ) can be interpreted as a macroscopic static residual obtained from 



the computation of a single-scale linear stage. In order to derive this term, the problem ( 13 1 must 
be solved independently on each substructure (the expression of the local macroscopic contributions 
{F_e)e<^'e is given in [[TO]]). The resolution of the macroscopic problem ([M]) leads to the global 



knowledge of Lagrange multiplier Wi , which is finally used as prescribed displacement to solve the 



substructure problems ( 13 1 



In order to perform the resolutions of ( 13 ) in substructure variables, finite element method is used. 



Since the constitutive law of the substructures is linear, the stiffness operator of each substructure can 
be factorized once at the beginning of the calculation and reused without updating throughout the 
analysis, which gives high numerical performance to the method. 

Algorithm [T] sums up the iterative procedure described in this section. 



Algorithm 1: The two-scale domain decomposition solver 



Construction of the stiffness operator of each substructure ; 

Computation of the macro homogenized operator of each substructure ; 

Global assembly of the macroscopic operator; 

Initialization sq G T; 

for n = 0, . . . , iV do 

Linear stage: computation of s„ G Ad ; 

Computation of the macro right-hand term Fe of each substructure ; 
Global assembly of the macroscopic right-hand term ; 



Resolution of the macroscopic problem ( 14 ) 
Resolution of the microscopic problems ( 13 1 
Local stage: computation of s„_|_i € T ; 



Resolution of the local problems on Interfaces {J^ ee'){e,e')i^'e'^'i 
Resolution of the local problems on Boundary interfaces {TEi)E£'E 
Computation of a global error indicator (distance separating Ad and F) 
end 
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2 Loss of scalability in the case of crack propagation 



Our first tests have outlined a numerical issue arising when simulating delamination propagation with 
the multiscale domain decomposition method. We observed a drop of the convergence rate which, as 



explained in subsection 2.1 can be interpreted as a loss of scalability and for which classical solution 
presented in subsection 2.2 are not efficient. Next section proposes a more realistic solution to this 
problem. 



2.1 Causes 



Test case 1 




perfect interfaces 
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macro forces 
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initial crack 



micro forces 



Figure 4: DCB-likc 2D test cases. Macroscopic and microscopic solution forces are represented in Case 
2. 



Figure |4] shows the results of a computation performed on DCB-like (double cantilever beam) 2D 
test cases. The constitutive law of the plies is linear and isotropic, while the interfaces are perfect 
bonds, except three of them which simulate an pre-existing crack (unilateral frictionless contact). In 
Test case 1, the prescribed traction leads to the computation of a homogeneous solution. The fields are 
perfectly represented in the macroscopic space, which results in a high convergence rate of the strategy 
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Test case 3 




micro forces cohesive interface 



Figure 5: DCB 2D test case (cohesive interface separating isotropic plies). The microscopic stress 
concentration around the tip of the crack is not transmitted to the far substructures by solving the 
macroscopic problem. 



(Figure [6]) . In the second case, a stress singularity appears at the tip of the crack. Unfortunatly 
the classical scale separation is not adapted to that shape of stress field: a significant part of the 
stress field is orthogonal to the classical macroscopic space though its zone of infiuence is not confined 
(see the micro forces in Figure |4|. Hence local resolutions (with communications between adjancent 
substructures) are used to transmit large wavelenght information resulting in a drastic drop of the 
convergence rate of the LaTIn solver (Figure |6]): in such a case the method is no more scalable. 

In the last of these 2D test cases, we introduce the debounding ability by replacing some of the per- 
fect bonds by cohesive interfaces (Figure [s]). The immediate effect is a further drop in the convergence 
rate (Figure [6]). Not only the stress singularity impairs the convergence (see micro forces distribution 
in Figure [5]) but also, in order the tip of the crack to propagate from one Gauss point to its neighbor, 
a sufficiently converged global equilibrium state is indeed required. Hence, the solution to debounding 
problems requires the computation of successive quasi-equilibrated states within each load increment, 
which penalizes the numerical efficiency of the strategy. 



2.2 Remedies 

Such a problem has already been encountered. In |5], the authors simulate a crack propagation 
within the LaTIn-based multiscale framework, using the fracture mechanics theory. The cracks are, 
in this case, cutting the interfaces of the domain decomposition strategy, which results in a drop 
in the convergence rate. Scalability is successfully restored by enriching the macroscopic basis with 
piecewise linear functions in order the discontinuous interface displacement fields to be represented in 
the macroscopic space. 

Following the same idea, we enrich the macroscopic basis in order to obtain a good representation 
of the stress concentration due to the cohesive crack in the macroscopic space. Though, in our case, 
the crack front position is not known in advance, and every potential position should be taken into 
account. Hence, no significant gain is obtained unless the macroscopic space is enriched locally around 
the tip of the crack (black square on Figure [5| by all the finite element interface functions (Figure [7]) . 
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10 20 30 40 

Iterations 



50 60 



Figure 6: Convergence curve of the DCB-like test cases. The error criterion is a normahzed distance 
separating Spaces and F |T2] 



This corresponds to clearing out the microspace and using maximal macrospace. 

The results obtained in terms of convergence rate are presented on Figure [8] During the considered 
increment, the propagation spreads along the distance covered by three Gauss points; the damage 
state of one Gauss point needs to be sufficiently well converged before the next one starts damaging 
(which justifies the swaying shape of the curves), once the damage state of Gauss points is correctly 
acquiered then the whole structure converges (fast decreasing part of the black curve) . The bump in 
the gray curve (classical macrospace) observed around Iteration 40 corresponds to the one observed 
around Iteration 8 of the black curve (full macrospace) . 
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Figure 7: Linear macroscopic basis and finite element interface functions in 2D. 



Though this enrichment is performed locally in small process zones, the extra computational costs 
are not affordable in 3D simulations. Indeed, computing the local homogenized operators (L^^)e£E 
requires to perform a number of resolution of Problems ( 13 1 equal to the number of interface degrees 
of freedom (explicit computation of Schur operators) . 

Moreover, one can easily show [lOj that such an enrichment, coupled with an adequate choice of 
the search direction E~ is equivalent to bond the subdomains with a linearized interface behavior 
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linear macro basis 
"full" enrichment 




Figure 8: Convergence curve of the DCB test case using the locaUy enriched macroscopic basis. 

equation Q at the hnear stages of the LaTIn solver. Hence, at this stage, a fully equilibrated solution 
is computed in the enriched zone. 

The remedy that we proposed in this section thus implies lots of computations which makes it only 
adapted to 2D problems. Indeed in 3D, the number of degrees of freedom within the damaging zone 
is too large to consider full condensation. Then in the next section we propose an alternative strategy, 
suited to 3D problems, to compute this local solution. 



3 Relocalization strategy in the vicinity of the crack 

In order to illustrate the relocalization technique, we consider the 3D version of the DCB problem 
previously studied (see Figure [9|). T he loading is applied in 10 increments, damage begins at the third 
one. As can be seem on Figure 10 ("No subresolution" data), as soon as damage appears the number 
of iterations per increment explodes (about 10 times more iterations). For such a problem, using a full 
macro basis would imply unrealistic extra computations. 



3.1 Principle of the sub-resolution strategy 

The alternative technique consists in extracting a part Qsub of Domain £7, where flgub is a source of 
localized nonlinearities in the structure. At each linear stage of the global LaTIn solver, the converged 
solution of the nonlinear subproblem in flgub is sought using the two-scale domain decomposition 
strategy described in Section [l] (see Figure ^) as together with Algorithm 

3.2 Local enhancement of the "linear" stage 

The new problem to solve at the enhanced "linear" stage now reads: Find s = {sE)EeE, where 
se = (We, Fe^ solution to: 

• linear equations in substructures variables in Domain Q 

— static and kinematic admissibility of the substructures, Equations ([T]) and ^ 

— linear constitutive law of the substructures. Equation ^ 
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Global LaTIn 
iterations 



interface which owns the 
maximum local damage rate 




extraction of a local subproblem 




prescribed robin boundary conditions 
(from the previous global computation) 



Figure 9: Illustration of the subresolution technique 



• interface bonding equations in ri\r2s„f, 

— search direction 

V TsE'/iE i n,ab or E' i ilsub), {s-s)e E" 

— equilibrium of the macroscopic interface forces 



F 



M 



V Tee'/{E i ilsub or E' i nsub), 
• interface bonding equations in ^sub'- interface linear and nonlinear behavior 

V TeE'I{E,E') € riLb, nEE'{WE,WE,.FE,FE,) - 

V TeJE e n,^b and Te n {dn^ ndQf)^ 0, REd{WE,FE) = 



(15) 

(16) 

(17) 
(18) 



Solving these equations requires to perform a classical linear stage on the whole structure (as 
described in Sect. [T]), and a nonlinear "exact" subresolution on i^sub- 

3.3 Nonlinear subresolutions 

On the boundary dflsub\dfl, the subproblem is bonded to the solution computed at the previous global 



linear stage. In order to enforce equations (15) and (16), macroscopic equilibrium and microscopic 
Robin boundary conditions are enforced: 



V TsE'/iE e n,ub and E' ^ r!,.,^), 



e W", I {F_E + k~W_E - (Eeg + k^Wsa)) -W!^* dT = {) 



(19) 
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Algorithm 2: The subresolution strategy algorithm 



Construction of the stiffness operator of each substructure ; 

Computation of the macro homogenized operator of each substructure ; 

Global assembly of the macroscopic operator; 

Initialization sq G F; 

for n = 0, . . . , iV do 

Enhanced linear stage: computation of s„ G Ad ; 
Classical linear stage on E ; 
Subresolutions : ; 

Locate the process zones requiring subresolution ; 
Assemble the macroscopic subproblem ; 
for j = 0, . . . ,m do 

Subproblem local stage ; 
Subproblem linear stage ; 
Computation of a local error indicator ; 
end 

Local stage: computation of s„_|_i G T ; 
Computation of a global error indicator 
end 



where F^eg ^^"^ ^eg respectively the force and displacement fields and W^^ obtained at the 
end global linear stage on {rEE'){Een,^t and E'^n.^t)- 

The set of equations in flsub (substructures admissibility and constitutive law in flsut, interface 



behavior on flsub\d^sub, boundary conditions on dVlsub H dfl, bonding equations (19) on dflsub\d^) 
is solved using the LaTIn two-scale procedure described in Section 1. Two particular technical points 
should be focused on: 



• The bonding equations p9| are handled at the local stage of this scheme in the same way the 
boundary equations ^ are dealt with in the classical two-scale resolution [T2]. However, these 
two equations, coupled with the search direction E+, lead to the resolution of two global systems 
on each interface ^EE'\{Een^^t and E'^n^^t)' ^'-'^ projections in the microscopic and macroscopic 
spaces are required. 

• At the linear stage of the subresolution scheme, the macroscopic equilibrium of Osub in enforced 
(so as to ensure the scalability of the subresolution strategy): 

V TEE'/iE G Qsub or E' G l}.„fc) , F%' + F%', ^ 

V TEE'/iE G n,ub and E' ^ Cisub), - Z|g = « (20) 

yvEjiEe andr^nar!/ 7^0), = o 

The resulting macroscopic subproblem reads: 



yw G w;„6 , 

, ^ f — M — M* . ^ f — M* 



rBn(an„,6\9n) (21) 



/ F^.W dV- Fe.W dV 

The construction of this problem is cheap as the homogenized operators L^^ are the same used for 
the global macroscopic problem (no enrichment of the macroscopic basis). Though, if mes(c}r2subn 
dilu = 0) , the macroscopic assembled operator has a Kernel (rigid body motions of the structure 
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extracted to perform the subresolution). Hence, Problem (21 ) must be solved using a generalized 
inverse. 



3.4 Results 



The example being treated is the DCB case represented Figure |9] In order to extract the subproblem 
automatically, we choose to perform the subresolution on a set substructures and interfaces in a box 
surrounding the cohesive interface with the highest damage rate at the end of the previous global 
LaTIn iteration (see Figure ( [q])). 

As it is shown on FigureflOl the resolution of a subproblem around the crack's tip leads to a 
convergence rate of the global LaTIn solver which is independent on the load increment of the analysis 
(i.e. independent on the area of the interface which becomes delaminated in one increment), which 
means that the numerical scalability is restored. In addition, the numerical complexity of the strategy 



is considerably reduced. In our case, the number of local inversions (resolution of problem (13l) is 
divided by two. 



250 



- 200 
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Figure 10: Subresolutions around the tip of the crack: results 



These results are obtained for a rather "large" subresolution box. The influence of the size of the 
subproblem on the numerical scalability is shown on Figure ( 11 ). One can clearly see that the number 



of global LaTIn iterations to convergence drops with an increasing size of the subresolution zone in 
an asymptotic manner (we observed that no significant gain is observed when increasing the size of 
the box defined in Case 3). Indeed, as soon as the global effects due to the stress concentration are 
enclosed in the subresolution zone, no further gain can be obtained by making use of this technique, 
thus the optimal size of the subresolution zone seems only to depend on local structural stiffnesses and 
dimensions (mechanical properties of the process zone) and not on the remaining of the structure. 



4 Discussion 

The global CPU time is not significantly reduced when implementing this strategy on parallel architec- 
tures. Indeed, using the initial allocation among the parallel processors addresses the subresolution to 
a very small number of processors. The first solution to tackle this difficulty is to perform independent 
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Time step 

Figure 11: Influence of the size of the subproblem 

subiterations systematically on all the processors [H . This strategy has successfully been assessed 
on various tests though its domain of efficiency is not well mastered yet. An other solution is to change 
the allocation of the substructures among the processors on the fly during the computation in order 
to take into account the level of nonlinearity of each set of substructures so that a balanced CPU load 
is reached. 

The second point to discuss about is the lack of generalization of this dedicated technique. Figure [12] 
clearly illustrates the difficulties that might arise when using the subresolution algorithm in the general 
case, for multiple crack fronts propagation may be involved (refer to [10] for details on the computation 
of such large debounding problems on parallel computers) . The front have an arbitrary shape, which 
raises the complex issue of the choice of the number and size of the relocalization zones and the problem 
of potential interactions between relocalization zones. In addition local instabilities might also appear 
during the relocalization computations, which is not accounted for at this stage of our developments. 




Figure 12: Damage map in the interfaces of a [0 ± 45 90]^ composite holed plate 

Hence, in the future, the relocalization strategy should be made more general and robust in order to 
improve the efficiency of the enhanced multiscale domain decomposition strategy for complex laminated 
structures. 
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5 Conclusion 

The accurate prediction of delamination in large process zones of laminated composite structures 
requires refined models of the material behavior. Such descriptions lead to the resolution of huge 
systems of equations. In order to compute the exact solution of such a refined model, we used a 
multiscale domain decomposition strategy based on a LaTIn iterative resolution algorithm. This 
method is particularly well-suited to solve problems where nonlinearities are introduced in surface 
joints bonding sets of 3D linear entities. 

We have shown that the classical scale separation was insufficient in the high gradient zones to 
provide numerical scalability. Using an enrichment approach, we proved that a microscopic resolution 
was required in the process zone at each iteration of the iterative solver. Following this observa- 
tion, we developed a subresolution procedure which preserved the numerical scalability of the crack 
propagation parallel calculation. This procedure still needs automation for complex structures (espe- 
cially concerning the shape and the number of relocalization box(es)) and regulation in case of local 
instabilities. 
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